Pair density waves and vortices in an elongated two-component Fermi gas 
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We study the vortex structures of a two-component Fermi gas experiencing a uniform effective 
magnetic field in an anisotropic trap that interpolates between quasi-one dimensional (ID) and 
quasi-two dimensional (2D). At a fixed chemical potential, reducing the anisotropy (or equivalently 
increasing the attractive interactions or increasing the magnetic field) leads to instabilities towards 
pair density waves, and vortex lattices. Reducing the chemical potential stabilizes the system. We 
calculate the phase diagram, and explore the density and pair density. The structures are similar 
to those predicted for superfluid Bose gases. We further calculate the paired fraction, showing how 
it depends on chemical potential and anisotropy. 
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Introduction — Quantized vortices play an essential 
role in understanding the behavior of type-II supercon- 
ductors and superfluids such as ^He. In cold gases, these 
vortices were the smoking gun for superfluidity [1]. Here 
we study how confinement influences the vortex struc- 
tures in a trapped gas of ultracold fermions. We use 
the microscopic Bogoliubov-de-Gennes (BdG) equations, 
and consider anisotropic traps that interpolate between 
quasi-one dimensional (ID) and quasi-two dimensional 
(2D). 

The behavior of topological defects in confined geome- 
tries can be quite rich. A good example is rotating bosons 
in anisotropic traps [2], where one sees multiple transi- 
tions in the structure of vortex lattices as the parame- 
ters are changed. Most intriguing, in the quasi- ID limit 
one sees a "roton" spectrum which softens as the ro- 
tation rate increases, signaling an instability to form a 
snake-like density wave. With recent experimental de- 
velopments [3], we expect these structures can soon be 
explored in Bose gases, and related studies will be under- 
taken in Fermi gases. In the Fermi gas, we find parallels 
to all of the predicted boson physics. The single particle 
instability which drives density waves in the Bose case 
becomes a collective instability for the fermions, and in- 
stead drives pair density waves [4] . For a range of param- 
eters we even find that the order parameter has the form 
predicted by Larkin and Ovchinnikov [5] for a polarized 
gas. 

In very different contexts, studies of vortices in con- 
fined geometries lead to a number of interesting and im- 
portant results such as "non-Hermitian" quantum me- 
chanical analogies [6], and the destruction of supercon- 
ductivity via phase slips [7]. Generically, reducing the 
dimensionality enhances fluctuations, leading to novel ef- 
fects. 

Driven partially by increased computer power and par- 
tially by interest in the BCS-BEC crossover, a number of 
research groups have recently produced Bogoliubov-de- 



Gennes (BdG) or density functional calculations of sin- 
gle vortices [13], and vortex lattices [14]. These have 
largely been 2D or three dimensional (3D) calculations, 
with translational symmetry along the magnetic field. 
The numerical challenges of these calculations come from 
the large basis set needed to describe the single particle 
states. By truncating to the lowest Landau level, one can 
greatly simplify the problem [15]. As we explain below 
this limit is experimentally relevant [16, 17]. 

Model — We start from the Hamiltonian of a spin bal- 
anced two-component Fermi gas, with a total number of 
particles N = J(^l^| + ^t^|)(ir and chemical potential 



/i. 



/C = 



E 



^tHo^. + Hint]dr-jlN, 



(1) 



where the single particle Hamiltonian Hq = {px — 
By)'^/2m + Py/2m + p1/2m + V(r), describes a neu- 
tral atom of mass m and momentum p experiencing 
a uniform effective magnetic field B in the z direction 
(Landau gauge), where the harmonic trap is V{r) = 
^m{ujyy'^ + cj^z^), and the inter-component interaction 

J^int = — ^^l^l^i^^, is attractive, with the coupling 
constant related to s-wave scattering lengths a^ via g = 
—ATih?as/m {g > 0) [10]. We do not treat the case where 
^ < 0, in which the physics is more involved [11]. This 
single particle Hamiltonian is readily engineered in cold 
atoms either by using two counter-propagating Raman 
beams with spatially dependent detuning [3] or rotating 
the gas in anisotropic traps where the rotation rate ap- 
proaches the weakest trapping frequency [12]. When ujz 
is large, this model can be tuned from quasi- ID to quasi- 
2D by changing ujy. 

(A) lowest Landau level — Following Sinha et.al 
[2], the single particle Hamiltonian is readily diagonal- 
ized, with eigenstates labeled by three quantum numbers 



K^n,n\ and energies given by 



(2) 



where the effective cyclotron frequency is ujc = 
j+cj^, the cyclotron frequency is uJc = B/m^ the 



characteristic energy of motion in the x direction is 
£ = hWy/Aujc^ and we have neglected the zero-point en- 
ergy. The dimensionless wave-number K = ^/2£k labels 
the momentum k along the x direction, where the ef- 
fective magnetic length is ^ = ^yh/mO^. The discrete 
quantum numbers n and n' corresponds to the number 
of nodes in the z and y directions. In the absence of 
confinement in the y direction, £ ^ 0, and we recover 
degenerate Landau levels. Hence, we refer to n as the 
Landau level index. If the interaction energy per particle 
(^int/^) and the characteristic "kinetic energy" {£K'^) 
are small compared to hwc and ficj^, one can truncate to 
the lowest eigenstates with n = n' = 0, which are of the 
form 



(t)K{r) 



1 i^^- 



yiridzL 






{y-VKr 



(3) 



where yK = V2uJcKI/2Cjc^ dz = \JhlvnijJz and L is the 
length in the x direction. 

The conditions allowing us to truncate to the lowest 
Landau level constrain the 3D density U'^b and magnetic 
field strength B. For example, the condition (Hint/^) <^ 
hwc ~ J^z requires usd ^ ^c/ 9 ^ ^z/ 9- The other 
condition, £K'^ <^ hujc^ requires B :^ muOy. While such 
fields are challenging to produce in cold atoms, they 
are not completely unreasonable. In a very recent ex- 
periment performed by I.Bloch's group [17], the den- 
sity is n^D ^ lO^^cm"^, and the cyclotron frequency is 
UJc ^ lOOkHz. Since this experiment involves coupled 
"wires", it is natural to use them for quasi- ID studies. 
Note, the magnetic field is "staggered" in that experi- 
ment, while we consider the uniform case. 

Letting ax annihilate the state in Eq. (3), one has an 
effective ID model, 

H/£ = Y,{K^ - fi)a^j,^aK. + /3 ^ fH<l)fiq), (4) 



K,( 



-l/8{2K-qy 



dq-KiCLK't^ the dimen- 



where f{q) = Y.k^ 
sionless chemical potential is /i = ji/£ and the effective 
interaction parameter is /3 = -|px(^)^^^(^)^- From 
the definition of /3, one sees that increasing the interac- 
tion strength g has the same effect as increasing the mag- 
netic field 5, increasing the ^-confinement uOz-, or reduc- 
ing the ?/- confinement ujy. In the following, we will inves- 
tigate the properties of the confined Fermi gas by study- 
ing Eq.(4). One can show that the interaction in Eq.(4) is 
equivalent to /^Xlg /U^)/(^) = f3 J drF^{r)F{r), where 



(B) Bogoliubov de Gennes approach — We intro- 
duce the pair field A^ = P{f{q))^ and its transform 
A(r) = P{F{r)). We neglect the fluctuation (/t(g) - 
Aq/P){f{q) - Aq/P) to reduce Eq.(4) to a bilinear form, 

H/£ = Y.{K^ - fi)a^^,aKa 

K,a 

+ ^(A*/(g) + Ajt(q)_|A,|V/3). (5) 



Given Ag, one can diagonalize H, and then impose 
self-consistency. For arbitrary Ag, this process is un- 
wieldy [14]. We here introduce two approximations which 
make the numerical calculations more efficient. First, 
we assume A^ is non- vanishing only when the central 
momentum of the paired fermions is q = nKo, where 
n = 0, ±1, ±2, .... The characteristic wave- number Kq is 
taken to be a variational parameter. This is equivalent to 
assuming A(r) is periodic in the x direction and treating 
the wavelength variationally. Second, we restrict our- 
selves to consider the symmetric pair field: A^ = A_g. 
This implies a spatially symmetric field A(r) = A(— r). 
Under these assumptions, the Hamiltonian is reduced to 



K,a n 

+ E(^Ni^o/(^^o) + ^Hifo/^(«^o)). (6) 



Since Eq.(6) will be calculated by taking the continuum 
limit ^^ -^ {'\/2L / ^irt) J dK [see Supplemental ma- 
terials)^ it is useful to introduce a positive parameter 
a = —V2LI3/4:7t£ to characterize the effective attrac- 
tive interaction. For small a, we find A\n\Ko 7^ for 
only a few values of n. We define ^ to be the number 
of nonzero A\n\Ko' The various phases can be distin- 
guished by looking at the pair density |(^^^|)p and/or 
the particle density (^1^^) (see Fig. 1(b)). The features 
are clearest in the pair density. If more than one AnKo 
is nonzero, we have either a pair density wave or vor- 
tices. For example, the case ^ = 3 (Aq ^ 0, A±Xq ^ 0), 
as illustrated in Fig. 1(b), corresponds to a pair density 
wave where |(^^^^)p has corrugations. The case ^ = 2 
(Aq = 0, A±Xo 7^ 0)7 consists of a single row of vor- 
tices. Larger (f, for example in Fig. 3, corresponds to a 
vortex lattice. The case ^ = 2 gives an order param- 
eter which can formally be identified with the Larkin- 
Ovchinnikov (LO) state [5] (see also [8]). Here, A^ is 
nonzero except when K = ±i^o- Defining an effective 
ID order parameter A^^{x) = ^^e*^^Ax, we have 
A^^{x) = 2Axo cosi^o^- Note that unlike the LO state, 
the physical order parameter A(r) = J^^ Ax^x(r), is 
not a simple cosine. Also note that unlike LO's model, 
here we assume both spin states have equal chemical po- 
tentials. Instead of being driven by the polarization, our 



instability towards a paired density wave is driven by the 
form of the effective ID interaction. 

When ^ = 1 (Aq 7^ 0), Eq.(6) can be analyzed ana- 
lytically {see Supplemental materials - A). One readily 
obtains the gap equation, 



(a) 



1 
a 



-K' 



2eK 



-dK, 



and the number equation, 
V2L 



N 



V2L r 

Ani J 



eo 



(1 - —)dK. 



eK 



ij) 



(8) 



/i. 



where ex = ^/^T~\A^^^^-^ and eo = K'^ 

Unlike the traditional case, the integrand in the RHS 
of Eq.(7) has a factor e~^ in the numerator, which dom- 
inates the behavior of the integrand for K ^ 1. If jj. ^ 1 
(meaning in physics units jj. ^ S)^ and Aq is sufficiently 
small, the integrand in Eq.(7) is bimodal. There is a gen- 
tle peak of height j- and width 1 centered at i^ = 0, and 



2^l 



a sharp peak of height 
atK 



-m/2 



and width 



|Ao|€ 



-m/2 



centered 



2|Ao| ^^^^ ^^^^^^^ Vm 

y7^. The power- law tails of this sharp peak give a 
contribution to the integral which scales as A^-^ log | Ao| 
as Aq ^ 0, where A is a constant. Solving Eq.(7) in this 
regime yields an extremely small order parameter. In 
this weak pairing limit, our numerics are unstable and 
the vortex lattices are better treated by expanding the 
energies in power of Aq [9] . 

Another instructive limit is /i < and N/L -^ 0, 
where the behavior is dominated by two-body physics. 
Eq.(7) then becomes the Schrodinger equation of a 
two-body problem in momentum space [18], i.e., a = 
2/ J e~^ /(^^ ~ fi)dK, where the two-body binding en- 
ergy T] is identified with twice the chemical potential. 

Phase diagram — We numerically minimize the en- 
ergy by studying Eq.(6) {see Supplemental materials - 
B). We find discrete jumps in ^ as a function of the di- 
mensionless attractive interaction a and the dimension- 
less chemical potential fi. The resulting phase diagram 
is shown in Fig. 1(a). The darkest red region (^ = 0) 
is the vacuum with no particles. Increasing a and/or /i 
brings one to a quasi- ID superfluid state. This state, 
characterized by ^ = 1, has no vortices and is transla- 
tional invariant in the x direction. The ^ = to (f = 1 
transition is continuous with Aq ^ and N/L ^ at 
the boundary. Further increasing a and/or /i leads to 
an instability towards a ^ = 3 state (the narrow yellow 
region). This state breaks translational symmetry. The 
transition is continuous, and the boundary can be found 
via a linear stability analysis of the ^ = 1 state {see Sup- 
plemental materials - C). At larger a and/or /i, there is 
a discontinuous transition to a state with ^ = 2. This 
sequence of instabilities closely mirrors what is found in 
calculations for Bose gases [2]. 
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FIG. 1: (color online) (a): The structure of phase diagram 
as a function of a and fi. The value of ^ (the number of 
nonzero A\n\Ko) is denoted in each region. The two black 
solid curves are the boundaries of two continuous transitions: 
^ = ^ ^ = 1 and ^ = 1 ^ ^ = 3. They show a fairly good 
agreement with numerics, (b): The structures of pair density 
I (^^^4.)!^ and density (^1^^) in the corresponding regions. 
The color key is shown in Fig. 3. 



Pair fraction — It is useful to put these results in the 
context of the BCS-BEC crossover. In 3D Fermi gases 
one thinks of the superfluid with /i < as being formed 
from tightly bound bosonic pairs, analogous to ^He. The 
superfluid with /i > is instead thought of within a BCS 
picture where diffuse pairs are formed by atoms at the 
Fermi surface. One can continuously tune between these 
two idealized limits by taking /i through zero: the size of 
the pairs varies continuously. Our approach to gaining 
insight into analogies with the 3D BCS-BEC crossover 
is to study the pair fraction P = 2A/'pair/iV [22], as in 
Fig. 2. While some of the qualitative features of the 3D 
crossover persist in our effective ID model, many of the 
details differ. 

To understand this figure, one must note that in a 
quasi- ID system the ratio of the interaction to the ki- 
netic energy is inverse proportional to the density, thus 
the strongly interacting regime can be reached by mak- 
ing the density small, or by making a large. The density 
increases monotonically with /i, but varies in a more com- 
plicated fashion with a. For small a and /i > we find 
dN/da < 0, while for large a and/or /i < we find 
dN/da > 0. At fixed a, the pair fraction decreases with 
/i (consistent with dN/dja > 0). 



1.0 



Low deasity 
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FIG. 2: (color online) The pair fraction P = 2Npair/N versus 
a with fi — —1,0, 1. The exponential small P for /i = 1 at 
a ^ is reminiscent of the BCS limit, and the large value of 
P for fi = —1 at a~ 1.5 is analogous to the BEG limit. The 
kink on each curve corresponds to the ^ = 3 ^ ^ = 2 phase 
transition. 




FIG. 3: (color online) The profile of density (left panel) and 
pair density (right panel) at a = 65, /i = 2, where the dimen- 
sionless coordinates are X = x/V2i, Y = y/y/2i. The color 
key is shown on the top. 



The top curve in Fig. 2, representing /i = — 1, starts at 
P = 1, roughly when a = 1.5. Such a large value of P is 
reminiscent of the BEG limit. The density vanishes here, 
then grows as a increases. For /i = — 1, the pair fraction 
decreases with a, except for a small kink, corresponding 
to the first order ^ = 3 ^ ^ = 2 phase transition. 

On the contrary, for /i = 1, P grows with a. As a ^ 0, 
P becomes exponentially small, as is predicted by the 
BCS theory. After a sharp rise, driven both by increasing 
a and decreasing N^ the pair fraction levels out. 

Each curve displays a kink, corresponding to the ^ = 
3 ^ ^ = 2 phase transition. As a increases to the region 
^ = 2, one row of vortices enters the elongated superfluid. 
This transition is accompanied by density modulations. 

To summarize we find that for /i > and small a 
the system behaves analogously to the BCS limit, while 
for /i < and a ^ |/i| the system behaves more like 
the BEG limit. The density vanishes if /i < and a < 
|/i|. For most of our parameter range, we observe physics 
analogous to the crossover regime. 

Vortex lattice — With increasing a, the number of 
Fourier components ^ increases, and the width in the y 
direction grows. We illustrate the large a limit in Fig. 3 by 
calculating the density and the pair density of the state 
with /i = 2,a = 65 and ^ = 7. Only "faint" vortices 
are seen in the density (left panel). Unpaired fermions 
fill the vortex cores leading to very poor contrast. On 
the contrary, one sees a clear stretched triangular lattice 
in the pair density (right panel). The lattice spacing 
is ~ 27T^/2i/Ko and the size of the vortex core is ~ i. 
Note the dimensionless wave-number Kq varies slightly 
with a but is of order 2. The vortex lattice is slightly 
deformed from a regular triangular lattice, but we expect 
this deformation to disappear in the quasi-2D limit {a -^ 
oo). 



Observation — Since the density depletion in the vor- 
tex core is highly suppressed, directly imaging the vor- 
tices through phase contrast or absorption imaging would 
be challenging. Coherent Bragg scattering of light may 
be a promising route for increasing the sensitivity of such 
optical probes [20]. One can also study the structures 
of pair density through photoassociation [21], where the 
paired state is transformed to a bound molecular state 
after illuminated with light. 

Summary — We have studied the two-component 
Fermi gases in elongated geometries. Truncating the 
BdG equations to the lowest Landau level, we investi- 
gate the vortex structures that emerge as the trap evolves 
from quasi-lD and quasi-2D. We calculate the phase di- 
agram and find instabilities towards pair density waves 
and vortex lattices. We explore the structures of den- 
sity and pair density, and calculate the pair fraction. We 
hope our results can soon be explored in experiment. 
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SUPPLEMENTAL MATERIALS 

A — Derivation of gap equation and number equation 

Here we analyze the special case where ^ = 1, corre- 
sponding to a ID model with translational invariance: 
Ag = unless (7 = 0. Under these circumstances, Eq.(6) 
simplifies to 

+ 53(A5/(0) + Ao/t(0))-|Ao|V/3- (9) 

K 

where we have introduced the dimensionless Hamiltonian 
Ho = H/S, and /(O) = ^j^ e~^^ a-KiaRt- 

Ho can be diagonalized in terms of non-interacting Bo- 
goliubov quasi-particle operators ^x, xk by the transfor- 
mation 






UK -Vk \ ( ^^ 



yielding the diagonalized Hamiltonian, 



K 



13 



where 



UK 



2€K 



-,VK 



ex - ep 
2eK 



(10) 



(11) 



(12) 



eK = ^Jel + \Ao\^e-^\eo = K'-^i (13) 

We introduce the dimensionless energy J^ = 
{2V27r£/ LS){GS\H\GS), where the ground state 
\GS) is annihilated by quasi-particle operators (,k-,Xk-, 



T= f{eo-eK)dK^\Aof/a, 



(14) 



where we have taken the continuum limit ^^ -^ 
{V2L/47Ti)JdK. 

Making (l/|Ao|)9J-'/9|Ao| = yields the gap equation 
(7). Letting N = -{V2L/47Tl)dT/dfi yields the number 
equation (8). These equations are further explored in the 
main text. 



B — Numerical approach 

We here describe our numerical approach to solving 
the BdG equations in the general case where ^ > 1. For- 
mally, Eq.(6) can be expressed in terms of non-interacting 
Bogoliubov quasi-particles by a canonical transformation 






E 



Vn' n, 






(15) 



where we have defined Kn = K — uKq^ and Un'n = 

Un'{Kn),Vn'n = Vn'{Kn). The matrix elements Un'n.Vn'n 

are governed by the following BdG equations, 



^Kr^ 



E 






F ,8 / 



(16) 



where e^^ is the dimensionless excitation energy of Bo- 
goliubov quasi-particles, and £n = K^ — /i, A^, = 

A|n|Ko6~^*^^^'''~^^°^ 5 ^^^ ^nn' IS the (5-function. In 
terms of the Bogoliubov operators the Hamiltonian is 
diagonal, 

r Ko/2 

n ^K^-Ko/2 
Ko/2 



K=-Ko/2 



.(17) 



The dimensionless ground state energy F = 
{2V27r£/LS){GS\H\GS) can be written as 

/ rKo/2 \ 

^ = E / (^n-exJ^i^+IAn^^jV^ . (18) 

^\J-Ko/2 ) 

For a given {/i, ifo, Ai^i^o}? we truncate Eq.(16), and use 
standard linear algebra packages to extract e^^. This 
effectively gives us J^ as a function of {/i, a, i^o, A|^|Xo}- 
This J^ is a variational upperbound on the true ground 
state energy. We fix {/i, a} and numerically minimize T ^ 
varying {iiTo, A|^|Xq}, using a quasi-Newton algorithm. 
We restrict the sum over n in Eq.(18) to — C < n < 
(". We find for the parameters studied, our results are 
independent of (" if (" > 6. 



C — Linear stability analysis 

Here we find the ^ = 1 to ^ = 3 phase boundary 
through a linear stability analysis. We take Aq > 0, 
and assume A^o = A_Xo = ^^ is small. We have chosen 



this factor of i, as the unstable direction will then yield 
real b. We will calculate D = d'^T/dS'^\s^o. For small a 
the curvature D is positive and the state with (5 = is 
stable. We find the instability by seeking the point with 
when D = 0. 

Within our ansatz for A\n\Ko^ the mean field Hamilto- 
nian is 



H/S = no ^iSA- 26^(3, 



(19) 



where 



A 



/t(Xo) + Pi-Ko) - fiKo) - fi-Ko). (20) 



Making use of the Hellmann-Feynman theorem, the 
second derivative of T can be expressed as 



05^ 



-i--^{GS\ppiKo)\GS) + - 



d5 



a 



(21) 



Setting D = d^T/d5^\s=o 
of instability is given by 



0, one finds that the points 



d 



i = pf^{GS\f\Ko)\GS). 



(22) 



Since the formal manipulations of perturbation theory 
are more transparent of finite temperature, it is conve- 
nient to rewrite Eq.(22) as 



— i = lim /3 



d Tr{e-^/^f\Ko)) 



r^o'^dS Tr{e-^/^) 



where T is a formal parameter. 

Substituting the results of Eq.(ll)-(12) to Eq.(23), we 
obtain 



"/ 



{uKUK-Ko^VKVK-Ko)^e 4(^0 ^X) 
^K-Kq + ^K 



■dK = 1.(24) 



This integral must be performed numerically, giving the 
second (right) black solid curve in Fig. 1(a). 



